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Aeroelastic Analysis of Modern Complex Wings 
Using ENSAERO and NASTRAN* 

Abstract 

A process is presented by which static aeroelastic analysis is performed using Euler 
flow equations in conjunction with an advanced structural analysis tool, NASTRAN. The 
process deals with the interfacing of two separate codes in the fields of computational 
fluid dynamics (CFD) and computational structural dynamics (CSD). The process is 
demonstrated successfully on an F/A-18 Stabilator (horizontal tail). 

Introduction 

With advanced subsonic transports entering into the transonic regime and fighter 
aircraft being limited by aeroelastic phenomena, it is becoming increasingly important to 
perform static and dynamic aeroelastic analysis using highly accurate fluid and structural 
models. Recently, work has been done using advanced structural finite element wing-box 
models coupled with Euler flow equations, but the structural finite element model is a 
simple one. 

There is a need to used highly detailed computation fluid dynamics (CFD) models 
coupled with highly detailed computational structural dynamics (CSD) models. Methods 
exist to solve each one independently, but not much work has been done coupling both 
highly accurate CSD and CFD models. 

This paper presents a process by which a highly detailed CSD model is coupled with 
a highly detailed CFD model. The analysis of the CSD model is done using NASTRAN, 
while analysis of the CFD model is done using Euler flow equations in NASTD (an in- 
house McDonnell Douglas Aerospace East CFD code). 

The aeroelastic process is broken down into the following steps: 

1) Get CFD solution (possibly rigid steady state) 

2) Calculate pressures at CFD grid points on the aerodynamic surface (Stabilator) 

3) Map pressures to forces on structural grid 

4) Obtain response of the structure 
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5) Map displacements to the CFD grid points on the aerodynamic surface 

6) Deform entire CFD grid 

7) Repeat steps until convergence criteria is met 

The above steps will sometimes be referred to later as one cycle. The above process 
will be demonstrated on an F/A-18 Stabilator (horizontal tail). 

The CFD grid topology of the F/A-18 Stabilator can be seen in Figure 1, while a more 
accurate view will be seen later. The grid was divided into two zones so that the CFD 
portion of the process could be run using PVM to reduce computational time. The zone 1 
grid contains the F/A-18 Stabilator, and the grid dimensions are 141 points around the 
perimeter of the airfoil, 69 points spanwise, and 61 points in the normal direction. Zone 2 
looks like a rectangular box whose dimensions are 3 1 points running streamwise, 69 
points running spanwise, and 121 points in the normal direction. The complete grid 
consists of more than 800,000 points, but the CFD code was run sequenced in each 
direction to reduce computational time. The k index for the tip of the Stabilator is 44, and 
the i index for the leading edge of the Stabilator is 71 . Also, zone 2 matches with zone 1 
as follows. The trailing edge of the Stabilator in zone 1 corresponds to i = 1 and i = 141, j 
= 1 , and k = 1 to 69 ( k = 44 is the Stabilator tip ). These points are the same as i = 1 , j = 
61 , k = 1 to 69 in zone 2. A detailed grid is shown for the upper surface of the F/A-18 
Stabilator in Figure 2. 

The finite element model of the F/A-18 Stabilator is shown in Figure 3. The model 
consists of approximately 2000 nodes, 12000 d.o.f. This is the model used by 
NASTRAN to form the stiffness matrix. In the next few paragraphs, we will describe 
some of the steps involved in the aeroelastic process in more detail. 

Step 1 involves getting a CFD solution. In this case, the rigid steady state solution 
was used. Next, the pressures were calculated using the CFD solution on the CFD grid. 

Step 3 maps the pressures on the CFD grid to forces on the structural grid. This 
involves a pre-processed mapping. The mapping consists of the following information. 
For each CFD point (i,k) on the surface of the Stabilator, the area on which the pressure 
acts and unit normal is calculated. Now, the magnitude and direction of the force due to 
unit pressure is known. The next step is to find a structural triangle that surrounds the 
CFD point (i,k). This can be difficult due to the irregular grids of some structural models. 

It is assumed that the structural grid is divided into an upper and lower surface 
structural grids with overlapping points possibly occurring at the leading and trailing 
edges and also at the tip. To find the structural triangle associated with the CFD grid 
point, the 20 closest structural nodes are found using the upper or lower surface structural 
grid depending on which surface the CFD point is located (Figure 4). Then all possible 
triangles using the 20 points are formed. Next, the triangles that do not contain the CFD 
point (i,k) as an interior point are eliminated (Figure 5). Then of the remaining triangles. 



the largest vertex distance is measured for each triangle, where vertex distance is the 
distance between the structural node and CFD point (i,k). Finally, the triangle with the 
smallest largest vertex distance is chosen. The triangularization of the structural domain 
to the CFD domain for the F/A-l 8 Stabilator can be seen in Figure 6. Now that the 
structural triangle is known, the area coordinates of the CFD point are used to distribute 
the force to the nodes of the structural triangle. So for each CFD point (i,k), the necessary 
weight factors and destination nodes are known in the pre-processing stage as well as the 
direction of the application of the loads. As a side note, the 20 closest points can be 
changed to 25 closest points depending on the density of the structural grid. 

Step 4 in the aeroelastic process is to solve the structural system of equations 
[K]{u g }={f}. The role of NASTRAN is a limited one, in that its only job is to produce 
the stiffness matrix in the pre-processing stage. In this example of the F/A-l 8 Stabilator, 
a reduced stiffness matrix is used. 

At this point of the process, the pressures on the CFD grid have been mapped to 
forces on the structural grid, and the structural system of equations have been solved, and 
now {u g }, the displacements are known. 

The next step is to deform the CFD grid according to the displacements on the 
structural model. To interpolate the structural node displacements to displacements on 
the grid points on the aerodynamic surface, a surface spline by Harder and Desmarais is 
used. The remaining grid is deformed according to a cosine spacing function. The CFD 
grid is deformed only in the z-direction (normal). The in-plane displacements are 
neglected. 

So to deform the surface grid, a surface spline technique was chosen. The surface 
spline system of equations become [A]{c}={u sp i}. The [A] matrix is dependent on the 
coordinates of the spline points. Therefore, once the spline points are chosen, the [A] 
matrix can be generated in the pre-processing stage, and {c} is the vector of coefficients 
of the surface spline, while {u spl } are the displacements of the spline points. 

The surface of the CFD grid which is deformed by the surface spline is shown in 
Figure 7. The red grid corresponds to zone 1 i = 1 to 141 , j = 1 , k = 1 to 69 ( k = 44 
being the tip). The blue grid corresponds to zone 2, i = 1 to 31 , j = 61 , k = 1 to 69. The 
spline points were chosen from the surface grid as follows. Eighty six points were chosen 
from the surface of the structural finite element model shown in Figure 8. Looking at 
Figure 7, it is seen that more points are needed to ensure smoothness of the surface grid. 
So, additional spline points were chosen at the boundary of the surface grid. Three rows 
of points running streamwise near the spanwise tip ( k = 67 to 69 ) in zone 1 and three 
rows of points near the streamwise tip ( i = 29 to 3 1) in zone 2 were chosen as spline 
points with zero displacement enforced as can be seen in Figure 9. Once the spline points 
are chosen, the [A] matrix is generated and stored in the pre-processing stage. 



Recalling, in the aeroelastic process, the displacements on the structural model have 
been calculated. The next step is to extract the spline point displacements {u sp ,} from the 
global displacements {u g }. Once {u sp |} are known, the system of equations [A] { c } = { u sp | } 
are solved, and the coefficients of the surface spline are known. 

The next step in the aeroelastic process is to deform the surface grid shown in Figure 
7. This is done using the surface spline equation which is a function of the x and y 
coordinates of the point. This is done easily since the coefficients of the surface spline 
are now known. 

Once the surface grid is deformed, the remaining grid (step 6) is deformed according 
to a cosine spacing function. It is enforced such that the outer boundaries corresponding 
to j =jmax in zone 1 and j = 1 and j = jmax in zone 2, do not move. So a spacing 
function that varies smoothly from 1 to 0 was is to deform the interior grid (Figure 10) 
which is itself a function of the normal index j. Once the surface grid is deformed, the 
remaining interior grid is deformed for each k = constant face. In order to prevent 
overlapping and crossed-sides of the CFD grid, a minimum spacing criteria is chosen, and 
the j = 26 face was found to meet the requirement. The points from j = 2 to 26 are 
deformed in each k = constant face the same amount as the j = 1 face. Then the spacing 
function is applied from j = 26 to jmax. This is done for zone 1 and zone 2. Zone 2 
indices are different, but the same deformation procedure is used. Figure 1 1 shows the k 
= constant face. The red grid is the grid falling within the j = 26 surface, while the blue 
grid is the grid exterior to the j = 26 surface. 

Finally, a convergence criteria is chosen. The L-2 norm of the CFD equations has to 
be converged and the Stabilator trailing edge tip displacement has to be converged. If 
not, then the process is repeated. Each time the process is repeated will be sometimes 
referred to as a cycle. Results are shown for the F/A-18 Stabilator at Mach 0.95, 1 degree 
angle of attack, at sea-level. 


Results 

Figure 12 shows the L-2 Norm of the zone 1 and zone 2 versus iteration. The rigid 
steady state solution was obtained after about 4000 iterations. The CFD grid was 
deformed and the CFD code was rerun starting with the previous step’s CFD solution. 
Each time the grid was deformed and rerun is here referred to as a cycle. About 9 cycles 
were performed in the F/A-18 Stabilator example. The L-2 norm shows good 
convergence of the CFD equations. Also, the integrated loads and moments are also 
shown in Figures 13-16. The spikes in the figures corresponds to the aeroelastic cycles 
where the grid is being deformed. There was not any regular interval between cycles. 
The loads can be seen converging towards the final flexible Stabilator solution. 

Figure 1 7 shows the displacements along the leading edge of the Stabilator. It can be 
seen that initially after the first cycle, the Stabilator deformed its greatest amount, and 
after cycle 2, it deformed its least amount. Figure 1 8 shows the convergence of the 
Stabilator tip airfoil section. This correlates well with Figure 17 also. The displacement 



of the trailing edge of the Stabilator tip is plotted versus cycle number in Figure 19. It 
shows the convergence of the displacements. 

The rigid and final flexible F/A-18 Stabilator is shown in Figure 20. The red solid 
figure is the rigid case, and the black grid lines are the deformed case. The deflections 
are scaled by a factor of 10 to help visualize the displacements. 

Figure 21 shows the pressure coefficient variation on the upper surface of the rigid 
and flexible Stabilator. Not a big variation can be seen by the plot, but the shock location 
on the upper surface moves forward a little in the flexible Stabilator case. Figure 22 
shows the Mach number variation on the upper surface of the rigid and flexible 
Stabilator. In this plot, it is easier to see that there is some change between the rigid and 
flexible cases, especially near the tip. There seems to be a discontinuity near the break of 
the trailing edge near the tip of the Stabilator in the flexible case. This is due to the fact 
that the CFD code was run sequenced, to reduce the computational time and effort. 

Conclusions 

The objective of this study was to create a process by which NASTRAN can be 
coupled with an advanced fluid model to perform aeroelastic analysis on highly detailed 
CFD and CSD models. The results show that this was achieved, but the process is 
somewhat problem dependent. The pressure to force mapping is a robust procedure since 
it requires only the structural node coordinates and CFD coordinates. It is independent of 
the type of elements used and does not require a rectangular domain. The displacement 
mapping is also good since it does not also require a rectangular domain and it’s only 
dependent on the coordinates of the spline points chosen. It does require a user to choose 
spline points in advance and thus becomes problem dependent. Overall, the process has 
been successfully demonstrated on the F/A-18 Stabilator. 



Edges Only View of CFD Zones 1 and 2 



Figure 1 



Upper Surface Undeformed CFD Stabilator Grid 
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Figure 3 



Mapping Pressures from CGD Grid to Forces on Structural Grid 



Figure 4 





Mapping Pressures from CFD Grid to Forces on Structural Grid 



Structural Node 




Mapping of CFD Points to Structural Triangles 




Surface Grid Used to Deform Interior Grid 
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Points used to Create Surface Spline 
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Cosine Spacing Function 
Multiplicative Factor vs. Normal Index 
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= Constant Face of CFD Grid (Zone 1) 



Figure 1 1 




F/A-18 Horizontal Stabilator 
Mach 0.95 Sea Level Static 
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F/A-18 Horizontal Stabilator 
Mach 0.95 Sea Level Static 
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F/A-18 Horizontal Stabilator 
Mach 0.95 Sea Level Static 
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F/A-18 Horizontal Stabilator 
Mach 0.95 Sea Level Static 
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F/A-18 Horizontal Stabilator 
Mach 0.95 Sea Level Static 
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Surface Spline Deflections 
F/A-18 Horizontal Stabilator 
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Wing Tip Displacement 
F/A-18 Stabilator 
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Deformed and Undeformed Upper Surface of Stabitator 
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Figure 21 





Figure 22 




